BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)
Ver com a Me:
- conferência dos recálculos, alpha FS parou de responder
- Olhar o código e resíduos, alguns já não parecem mais ok
- marcar na planilha os bichos que mudaram
- cálculo beta RC dentro do wrangling ou separado? onde jogo as funções (em script separado e usar source)
- sugestão de como plotar quando tem nulo (linha tracejada mais clara e sem SE?)
- Plot também vai pro Git? script separado?
Model formulas
For each of our response variables we run three candidate models for each scale:
Model 0 <- ~1, absence of effect
Model 1 <- percentage of forest cover at 3 and 5km radius
Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)
Loading data and gathering into one dataset| Landscape | gamma_FS | gamma_NFS | alpha_FS | alpha_NFS | betad_FS | betad_NFS | ab.land_FS | ab.land_NFS | brc.fsNA | brc.nfsNA | brc.abund.fsNA | brc.abund.nfsNA | fc_3km | fc_5km |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 148 | 11 | 9 | 3.75 | 2.25 | 7.25 | 6.75 | 66 | 40 | -0.2928571 | 0.1826667 | -0.0014658 | 0.4252252 | 38.52523 | 43.86337 |
| 215 | 19 | 7 | 5.62 | 2.00 | 13.38 | 5.00 | 174 | 51 | 0.1057143 | -0.3123810 | 0.2759545 | 0.3903427 | 30.15485 | 25.23160 |
| 263 | 16 | 13 | 6.00 | 4.50 | 10.00 | 8.50 | 393 | 214 | -0.4919048 | -0.0323810 | 0.6333476 | 0.3923447 | 10.10641 | 13.41180 |
| 266 | 11 | 8 | 3.25 | 1.38 | 7.75 | 6.62 | 70 | 12 | -0.2257143 | 0.1338095 | 0.1087754 | 0.1472425 | 30.87525 | 31.16090 |
| 282 | 17 | 11 | 6.50 | 4.12 | 10.50 | 6.88 | 219 | 212 | -0.3275000 | -0.1152381 | 0.2126412 | 0.7177654 | 15.26270 | 16.35308 |
| 291 | 9 | 3 | 3.71 | 1.71 | 5.29 | 1.29 | 150 | 49 | -0.3866667 | -0.5904762 | 0.1027218 | -0.4251871 | 45.74294 | 45.13213 |
| 317 | 12 | 4 | 5.00 | 1.38 | 7.00 | 2.62 | 152 | 56 | -0.3853571 | -0.5321429 | 0.1819677 | -0.4317889 | 48.75805 | 47.85138 |
| 329 | 21 | 14 | 8.62 | 5.75 | 12.38 | 8.25 | 369 | 262 | -0.0871429 | 0.2217857 | 0.9169527 | 0.4891320 | 14.73791 | 13.62588 |
| 333 | 14 | 15 | 4.50 | 4.62 | 9.50 | 10.38 | 158 | 219 | -0.1275000 | 0.1932143 | 0.6788932 | 0.6190476 | 13.71403 | 14.50960 |
| 335 | 15 | 8 | 5.75 | 2.38 | 9.25 | 5.62 | 140 | 113 | -0.2564286 | -0.3366667 | 0.3479193 | -0.9437533 | 23.87913 | 19.21670 |
| 359 | 16 | 17 | 6.00 | 5.50 | 10.00 | 11.50 | 149 | 133 | -0.2253571 | 0.3128571 | 0.7704490 | 0.8121693 | 18.10208 | 17.62431 |
| 399 | 20 | 14 | 7.88 | 5.12 | 12.12 | 8.88 | 432 | 169 | -0.1223810 | 0.2703571 | 0.9618190 | 0.6938009 | 29.73616 | 24.97851 |
RESULTS
We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).
FOREST SPECIALISTS
Gamma diversity
Modelling
g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g5k -27.9 64.8 4.7 0.0 3 0.731
## g3k -29.6 68.2 3.0 3.4 3 0.136
## g5kq -27.8 69.4 4.8 4.6 4 0.074
## g0 -32.6 70.6 0.0 5.8 2 0.041
## g3kq -29.2 72.2 3.4 7.3 4 0.019
Summary top model
summary(g5k)##
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.23181 -0.12743 -0.03948 0.15159 0.27856
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.68384 1.96964 10.50 1.01e-06 ***
## fc_5km -0.21451 0.05782 -3.71 0.00404 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.0352002)
##
## Null deviance: 0.74994 on 11 degrees of freedom
## Residual deviance: 0.34328 on 10 degrees of freedom
## AIC: 61.839
##
## Number of Fisher Scoring iterations: 3
Inspect residuals:
sim <- simulateResiduals(g5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1487, p-value = 0.632
## alternative hypothesis: two.sided
hnp::hnp(g5k)## Gamma model
boot::glm.diag.plots(g5k) ## for glm poisson or neg binomPrediction:
Alpha diversity
Modelling
a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a5k -19.3 47.6 2.8 0.0 3 0.523
## a3k -20.3 49.5 1.9 1.9 3 0.199
## a0 -22.1 49.6 0.0 2.0 2 0.195
## a5kq -19.0 51.8 3.1 4.2 4 0.064
## a3kq -20.3 54.2 1.9 6.6 4 0.019
Summary top model
summary(a5k)##
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.43224 -0.11774 -0.03674 0.08987 0.35754
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.42155 0.94502 7.853 1.38e-05 ***
## fc_5km -0.07200 0.02815 -2.558 0.0285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.06034959)
##
## Null deviance: 0.97549 on 11 degrees of freedom
## Residual deviance: 0.61269 on 10 degrees of freedom
## AIC: 44.608
##
## Number of Fisher Scoring iterations: 4
Inspect residuals:
sim <- simulateResiduals(a5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1043, p-value = 0.656
## alternative hypothesis: two.sided
hnp::hnp(a5k)## Gamma model
boot::glm.diag.plots(a5k)Prediction:
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b5k -22.1 53.3 5.1 0.0 3 0.690
## b5kq -21.4 56.6 5.8 3.3 4 0.133
## b3k -23.9 56.8 3.3 3.6 3 0.116
## b3kq -22.7 59.2 4.5 5.9 4 0.036
## b0 -27.2 59.8 0.0 6.6 2 0.026
Summary top model
summary(b5k)##
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.24293 -0.12664 -0.05680 0.08772 0.34333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.29031 1.24634 10.663 8.8e-07 ***
## fc_5km -0.14349 0.03623 -3.961 0.00268 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03508083)
##
## Null deviance: 0.76409 on 11 degrees of freedom
## Residual deviance: 0.32769 on 10 degrees of freedom
## AIC: 50.257
##
## Number of Fisher Scoring iterations: 4
Inspect residuals:
sim <- simulateResiduals(b5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1585, p-value = 0.608
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b5k) ## for glm Gamma or neg binomPrediction
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc0 5.4 -5.4 0.0 0.0 2 0.475
## rc3kq 8.7 -3.7 3.3 1.7 4 0.203
## rc5k 5.9 -2.7 0.5 2.7 3 0.124
## rc5kq 8.1 -2.6 2.8 2.9 4 0.113
## rc3k 5.5 -2.0 0.1 3.4 3 0.085
Summary of 2nd best model, as the top is of absence of effect (~1)
summary(rc3kq)##
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.13500 -0.09118 -0.03185 0.07770 0.23076
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6988850 0.2220491 -3.147 0.0118 *
## fc_3km 0.0413030 0.0173151 2.385 0.0409 *
## I(fc_3km^2) -0.0007386 0.0002921 -2.529 0.0323 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.135 on 9 degrees of freedom
## Multiple R-squared: 0.4268, Adjusted R-squared: 0.2995
## F-statistic: 3.351 on 2 and 9 DF, p-value: 0.08171
Inspect residuals:
sim <- simulateResiduals(rc3kq)
plot(sim)## qu = 0.25, log(sigma) = -3.759907 : outer Newton did not converge fully.
testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)## Gaussian model (lm object)
Prediction:
Total abundance
Modelling
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -70.5 150.1 2.1 0.0 3 0.399
## gab0 -72.6 150.5 0.0 0.5 2 0.318
## gab3k -71.2 151.3 1.4 1.2 3 0.217
## gab5kq -70.4 154.5 2.2 4.5 4 0.043
## gab3kq -71.0 155.8 1.6 5.7 4 0.023
Summary 2nd best model, as top model was the one of absence of effect (~1):
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.537178159,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6496 -0.9873 -0.2873 0.6438 1.8333
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.93023 0.31931 18.572 <2e-16 ***
## fc_5km -0.02485 0.01110 -2.239 0.0251 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.5372) family taken to be 1)
##
## Null deviance: 17.268 on 11 degrees of freedom
## Residual deviance: 12.411 on 10 degrees of freedom
## AIC: 147.08
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.54
## Std. Err.: 1.83
##
## 2 x log-likelihood: -141.082
Inspect residuals:
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1614, p-value = 0.56
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS",
"FC at 5km",
round(gab.fs.aic$AICc[1],2),
round(gab.fs.aic$dAICc[1],1),
gab.fs.aic$df[1],
round(gab.fs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "~ 1 (NULL)",
round(gab.fs.aic$AICc[2],2),
round(gab.fs.aic$dAICc[2],1),
gab.fs.aic$df[2],
round(gab.fs.aic$weight[2],2),
"-"),
c(" ", "FC at 3km",
round(gab.fs.aic$AICc[3],2),
round(gab.fs.aic$dAICc[3],1),
gab.fs.aic$df[3],
round(gab.fs.aic$weight[3],2),
round(gab3k.r2$R2_Nagelkerke,2)))Prediction:
bRaup-Crick (abundance-based)
Modelling:
### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc5k 0.2 8.6 3.8 0.0 3 0.593
## rc3k -0.7 10.4 2.9 1.8 3 0.237
## rc0 -3.6 12.5 0.0 3.9 2 0.085
## rc5kq 0.3 13.1 3.9 4.5 4 0.062
## rc3kq -0.7 15.2 2.9 6.5 4 0.022
Summary top model:
summary(rc5k)##
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39458 -0.18080 -0.00713 0.15192 0.50954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.900965 0.175354 5.138 0.000439 ***
## fc_5km -0.017963 0.006073 -2.958 0.014339 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2607 on 10 degrees of freedom
## Multiple R-squared: 0.4666, Adjusted R-squared: 0.4133
## F-statistic: 8.749 on 1 and 10 DF, p-value: 0.01434
Inspect residuals:
sim <- simulateResiduals(rc5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)## Gaussian model (lm object)
Prediction
NON-FOREST SPECIALISTS
Gamma diversity
Modelling:
### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.nfs.aic <- AICctab(g0,
g5k, g5kq,
g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g3k -27.3 63.5 7.6 0.0 3 0.7911
## g5k -29.2 67.5 5.7 4.0 3 0.1095
## g3kq -27.1 68.0 7.8 4.4 4 0.0864
## g5kq -29.2 72.2 5.7 8.6 4 0.0105
## g0 -34.9 75.2 0.0 11.6 2 0.0024
Summary top model:
summary(g3k)##
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.38392 -0.22774 -0.05573 0.15121 0.43058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.37185 2.13069 8.622 6.07e-06 ***
## fc_3km -0.30280 0.05227 -5.793 0.000175 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.08050867)
##
## Null deviance: 2.71637 on 11 degrees of freedom
## Residual deviance: 0.78116 on 10 degrees of freedom
## AIC: 60.54
##
## Number of Fisher Scoring iterations: 4
Inspect residuals:
sim <- simulateResiduals(g3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90839, p-value = 0.96
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(g3k) ## for glm poisson or neg binomPrediction:
Alpha diversity
Modelling:
### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a3k -16.7 42.3 5.7 0.0 3 0.490
## a5k -16.9 42.9 5.4 0.6 3 0.371
## a5kq -16.2 46.2 6.1 3.9 4 0.071
## a3kq -16.5 46.6 5.9 4.3 4 0.057
## a0 -22.3 50.0 0.0 7.7 2 0.011
Summary top model:
summary(a3k)##
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.68238 -0.15994 0.00108 0.11140 0.55352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.82215 0.83950 6.935 4.02e-05 ***
## fc_3km -0.09207 0.02113 -4.357 0.00143 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1163241)
##
## Null deviance: 3.1391 on 11 degrees of freedom
## Residual deviance: 1.2537 on 10 degrees of freedom
## AIC: 39.336
##
## Number of Fisher Scoring iterations: 6
Inspect residuals:
sim <- simulateResiduals(a3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.93258, p-value = 0.928
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(a3k) ## for glm poisson or neg binomPrediction:
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b3k -25.0 58.9 5.9 0.0 3 0.753
## b5k -26.8 62.7 4.0 3.7 3 0.116
## b3kq -24.6 62.9 6.3 3.9 4 0.105
## b5kq -26.6 67.0 4.2 8.1 4 0.013
## b0 -30.9 67.1 0.0 8.2 2 0.013
Summary top model:
summary(b3k)##
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.71147 -0.22545 -0.03265 0.18829 0.45483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.59989 1.68741 7.467 2.14e-05 ***
## fc_3km -0.21226 0.04078 -5.205 0.000399 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1116167)
##
## Null deviance: 3.2119 on 11 degrees of freedom
## Residual deviance: 1.2312 on 10 degrees of freedom
## AIC: 55.933
##
## Number of Fisher Scoring iterations: 4
Inspect residuals:
sim <- simulateResiduals(b3k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65877, p-value = 0.68
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b3k) ## for glm poisson or neg binomPrediction:
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -0.6 10.2 2.3 0.0 3 0.390
## rc0 -2.9 11.1 0.0 0.9 2 0.250
## rc5k -1.2 11.4 1.7 1.2 3 0.219
## rc3kq 0.5 12.7 3.4 2.5 4 0.114
## rc5kq -0.9 15.5 2.0 5.3 4 0.027
Summary of 2nd best model, 1st wqs of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32488 -0.21626 -0.05457 0.24366 0.39972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.322872 0.191341 1.687 0.1224
## fc_3km -0.014015 0.006519 -2.150 0.0571 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2787 on 10 degrees of freedom
## Multiple R-squared: 0.3161, Adjusted R-squared: 0.2477
## F-statistic: 4.622 on 1 and 10 DF, p-value: 0.05707
Inspect residuals:
sim <- simulateResiduals(rc3k)
plot(sim)## qu = 0.75, log(sigma) = -3.658389 : outer Newton did not converge fully.
testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Prediction:
Total abundance
Modelling:
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -63.9 136.9 5.3 0.0 3 0.404
## gab3k -64.3 137.6 4.9 0.7 3 0.284
## gab5kq -62.1 137.8 7.2 0.9 4 0.252
## gab3kq -63.7 141.2 5.5 4.3 4 0.047
## gab0 -69.2 143.8 0.0 6.9 2 0.013
Summary top model:
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.164123767,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8516 -0.4183 0.2175 0.4846 0.8634
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.93704 0.33661 17.638 < 2e-16 ***
## fc_5km -0.04838 0.01182 -4.094 4.24e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.1641) family taken to be 1)
##
## Null deviance: 29.266 on 11 degrees of freedom
## Residual deviance: 12.765 on 10 degrees of freedom
## AIC: 133.88
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.16
## Std. Err.: 1.77
##
## 2 x log-likelihood: -127.88
Inspect residuals:
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.45825, p-value = 0.288
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS",
"FC at 5km",
round(gab.nfs.aic$AICc[1],2),
round(gab.nfs.aic$dAICc[1],1),
gab.nfs.aic$df[1],
round(gab.nfs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "FC at 3km",
round(gab.nfs.aic$AICc[2],2),
round(gab.nfs.aic$dAICc[2],1),
gab.nfs.aic$df[2],
round(gab.nfs.aic$weight[2],2),
round(gab3k.r2$R2_Nagelkerke,2)),
c(" ", "Non-linear FC at 5km",
round(gab.nfs.aic$AICc[3],2),
round(gab.nfs.aic$dAICc[3],1),
gab.nfs.aic$df[3],
round(gab.nfs.aic$weight[3],2),
round(gab5kq.r2$R2_Nagelkerke,2)
))Prediction:
bRaup-Crick (abundance-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -7.4 23.8 2.0 0.0 3 0.381
## rc0 -9.4 24.1 0.0 0.2 2 0.337
## rc5k -8.0 24.9 1.4 1.1 3 0.219
## rc3kq -7.3 28.3 2.1 4.5 4 0.041
## rc5kq -7.9 29.5 1.5 5.7 4 0.022
Summary 2nd best model, 1st is of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24646 -0.18490 0.04466 0.26676 0.52334
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.84188 0.33775 2.493 0.0318 *
## fc_3km -0.02258 0.01151 -1.962 0.0781 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4919 on 10 degrees of freedom
## Multiple R-squared: 0.278, Adjusted R-squared: 0.2058
## F-statistic: 3.851 on 1 and 10 DF, p-value: 0.07814
Inspect residuals:
sim <- simulateResiduals(rc3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Prediction:
AICc-model selection table
| Diversity | Model | AICc | dAICc | df | weight | r2 |
|---|---|---|---|---|---|---|
| Gamma FS | FC at 5km | 64.84 | 0 | 3 | 0.73 | 0.55 |
| Alpha FS | FC at 5km | 47.61 | 0 | 3 | 0.52 | 0.38 |
| FC at 3km | 49.54 | 1.9 | 3 | 0.2 | 0.27 | |
| ~ 1 (NULL) | 49.58 | 2 | 2 | 0.19 |
|
|
| Beta FS | FC at 5km | 53.26 | 0 | 3 | 0.69 | 0.58 |
| RCbeta (P/A-based) FS | ~ 1 (NULL) | -5.44 | 0 | 2 | 0.47 |
|
| Non-linear FC at 3km | -3.74 | 1.7 | 4 | 0.2 | 0.3 | |
| FC at 5km | -2.75 | 2.7 | 3 | 0.12 | -0.01 | |
| Total abundance FS | FC at 5km | 150.08 | 0 | 3 | 0.4 | 0.44 |
| ~ 1 (NULL) | 150.54 | 0.5 | 2 | 0.32 |
|
|
| FC at 3km | 151.31 | 1.2 | 3 | 0.22 | 0.32 | |
| RCbeta (abund-based) FS | FC at 5km | 8.61 | 0 | 3 | 0.59 | 0.41 |
| FC at 3km | 10.44 | 1.8 | 3 | 0.24 | 0.32 | |
| Gamma NFS | FC at 3km | 63.54 | 0 | 3 | 0.79 | 0.74 |
| Alpha NFS | FC at 3km | 42.34 | 0 | 3 | 0.49 | 0.63 |
| FC at 5km | 42.89 | 0.6 | 3 | 0.37 | 0.61 | |
| Beta NFS | FC at 3km | 58.93 | 0 | 3 | 0.75 | 0.65 |
| RCbeta (P/A-based) NFS | FC at 3km | 10.2 | 0 | 3 | 0.39 | 0.17 |
| ~ 1 (NULL) | 11.09 | 0.9 | 2 | 0.25 |
|
|
| FC at 5km | 11.35 | 1.2 | 3 | 0.22 | 0.17 | |
| Total abundance NFS | FC at 5km | 136.88 | 0 | 3 | 0.4 | 0.82 |
| FC at 3km | 137.58 | 0.7 | 3 | 0.28 | 0.79 | |
| Non-linear FC at 5km | 137.82 | 0.9 | 4 | 0.25 | 0.93 | |
| RCbeta (abund-based) NFS | FC at 3km | 23.84 | 0 | 3 | 0.38 | 0.21 |
| ~ 1 (NULL) | 24.08 | 0.2 | 2 | 0.34 |
|
|
| FC at 5km | 24.94 | 1.1 | 3 | 0.22 | 0.13 |
FIGURES
Gamma hab
### Gamma ####
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g5k, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g3k, type="response", newdata=new2, se.fit=T)
new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96
new2$ic.low <- preds2$fit - preds2$se.fit*1.96
## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
xlab(" ") + ylab("Gamma") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=20.5, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.66"))) +
annotate("text", y=19, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.74"))) + labs(tag="a"))Alpha hab
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a5k, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
## prediction a1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a3k, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
xlab("") + ylab("Alpha") +
scale_color_manual(breaks=c("FS", "NFS"),
values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=14,face="bold"),
legend.position = 'bottom') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=13, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0383")))+
annotate("text", y=11, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.63"))) +
labs(tag="b"))Beta hab
### Beta add ####
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b5k, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
## prediction b1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b3k, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
geom_line(data= new2, aes(x= fc_3km, y= pred),
col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
xlab(" ") + ylab("Beta") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=18, x=40,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.58")))+
annotate("text", y=16, x=38, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.65"))) +
labs(tag="c"))## plot here not to need posterior figure editing app
g.hab + a.hab + b.hab#ggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)Total abundance FS & NFS
## FS
gab1 <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab1, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
gab1 <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
## prediction gab2
new2 <- new
preds <- predict(gab1, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) +
geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
# geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
scale_x_continuous(limits= c(10,50)) +
xlab(" ") + ylab("Total abundance") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
# scale_x_reverse() +
annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.44"))) +
annotate("text", y=270, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.82"))) +
labs(tag="d"))
### Beta RC FS & NFS
## FS
rc2 <- lm(brc.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc2, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
rc1 <- lm(brc.nfsNA ~ fc_3km, data=data_div)
## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) +
geom_line(data= new,aes(x=fc_3km, y=pred),
col= "#0072B2", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
legend.position = 'none') + ylim(-1,1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
label=expression(paste(~ R^2 ~ "= 0.82"))) +
annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.44"))) +
# annotate("text", y=0.9, x=43, colour= "#D33E00",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
# annotate("text", y=0.8, x=43, colour= "#0072B2",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
labs(tag="e"))
Beta RC abund FS & NFS
rc2 <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
rc1 <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(-1,1.1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=45, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.41"))) +
annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
label=expression(paste( ~ R^2 ~ "= 0.21"))) +
labs(tag="f"))Final plot
(plot.hab <- (g.hab + a.hab + b.hab) /
(gab.hab + rc.hab + rc.ab.hab))#ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)